
#The following script replicated analysis in Soroka, Fournier and Nir, "Cross-national evidence of a negativity bias in psychophysiological reactions to news," in PNAS.  Analysis are conducted on several datasets, as follows:
  
#datafile: workingdata1.Rdata
#dataframe: rs.hrv
#description: This is the respondent-stimulus dataset used for HRV analysis in Tables S2, S3 and S5, as well as Fig 1.
#variables:
#resp: respondent number
#country: respondent country
#stimulus: story identifier
#stimtype: stimulus type
#rmssd: RMSSD
#negativity: coder-attributed negativity score
#order: story order
#wt: country-level weights
#subtitles: =1 where respondent used subtitles

#datafile: workingdata2.Rdata
#dataframe: rs.gsl
#description: This is the respondent-stimulus dataset used for GSL analysis in Tables S2, S3 and S5.
#variables:
#resp: respondent number
#country: respondent country
#stimulus: story identifier
#stimtype: stimulus type
#gslmc: normalized GSL
#negativity: coder-attributed negativity score
#order: story order
#wt: country-level weights
#subtitles: =1 where respondent used subtitles

#datafile: workingdata3.Rdata
#dataframe: p.gsl
#description: This is the time-series panel dataset used for GSL analysis in Tables S4.
#variables:
#resp: respondent number
#country: respondent country
#country2: respondent country, distinguishing between country subsamples
#story: news story number
#timesec: seconds, over entire experiments
#timesecStory: seconds, counting by story
#timesecStorylog: logged timesecStory
#order: story order
#gslmc: normalized GSL
#negativity: coder-attributed negativity score

#datafile: workingdata4.Rdata
#dataframe: sv
#description: This is the survey (by-respondent) dataset used for descriptives in the Appendix, and for Figures 3 and 4.
#variables:
#country2: country, distinguishing between country subsamples
#respnum: respondent number
#gender: respondent gender
#country: country
#age: respondent age
#story11.rating3: negativity rating for story on Chimbolte
#story12.rating3: negativity rating for story on May Day
#story13.rating3: negativity rating for story on Niger
#story14.rating3: negativity rating for story on Sir Lanka
#story15.rating3: negativity rating for story on Gorillas
#story16.rating3: negativity rating for story on Folding Car
#story17.rating3: negativity rating for story on Young Director
#story18.rating3: negativity rating for story on Cure Cancer
#negcoef.rmssd: respondent negativity coefficient, for RMSSD
#negcoef.gsl: respondent negativity coefficient, for GSL
#dropgsl.new: (revised) variable indicating flawed GSL data
#dropbmp: variable indicating flawed BMP data
#hofstede: country-level measure, Hofstede's uncertainty avoidance
#gdppc: country-level measure, GDPPC
#media: country-level measure, Media Freedom

#Note that due to missing data (errors with physioligical sensors), the number of cases varies slightly across files. Country-level weights are thus assigned separately for each dataset.

#Note also that some analyses and all graphics are done in R, but the panel estimations are done in STATA. We do so through R, using the RStata package, below. Note all tests of statistical signficance are produced below, although they should be straightforward to produce with the information that follows.  All the critical regression models, as well as the ridgeplots, are generated below.

#####
#setup
library(stargazer) ; library(lmtest)
library(DataCombine) ; library(memisc)
library(ggridges) ; library(ggplot2)
library(RStata)
options("RStata.StataPath" = '/Applications/Stata/StataMP.app/Contents/MacOS/stata-mp')
options("RStata.StataVersion" = 15)

#####
#Loading datasets
setwd("/Users/stuartsoroka/Dropbox/A-drive/_Research/_Physiology/_Comparative/_papers/_pnas")
load("workingdata1.Rdata")
load("workingdata2.Rdata")
load("workingdata3.Rdata")
load("workingdata4.Rdata")

#####
#HRV results in Table S2 and Table S3

stata_src <- "
  xtset resp order
  quietly xtreg rmssd negativity order , fe 
  estimates store A
  quietly xtreg rmssd negativity order [weight=wt], fe 
  estimates store B
  esttab A B, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w r2_b r2_o N N_g)
  quietly xtreg rmssd negativity order , fe 
  margins, at(negativity=(-2 0 2)) post
  "
stata(stata_src,data.in=rs.hrv)

#####
#GSL results in Table S2 and S3

stata_src <- "
  xtset resp order
  quietly xtreg gslmc negativity order , fe 
  estimates store A
  quietly xtreg gslmc negativity order [weight=wt], fe 
  estimates store B
  esttab A B, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w r2_b r2_o N N_g)
  "
stata(stata_src,data.in=rs.gsl)

#####
#GSL results in Table S4 and Figure 2

stata_src <- "
  xtset resp timesec
  gen lgslmc = L.gslmc
  quietly xtreg D.gslmc c.negativity##c.timesecStorylog order lgslmc, fe 
  estimates store D
  esttab D , cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w r2_b r2_o N N_g)  
  margins, at(negativity=(-2 0 2) timesecStorylog=(2.99))
  "
stata(stata_src,data.in=p.gsl)

#####
#HRV results in Figure 3

c <- sv[sv$dropbmp<2,c("negcoef.rmssd","country")]
c <- na.omit(c)
c$countryt <- NA
c$countryt[c$country=="GH"] <- "Ghana"
c$countryt[c$country=="RU"] <- "Russia"
c$countryt[c$country=="BR"] <- "Brazil*"
c$countryt[c$country=="CH"] <- "Chile"
c$countryt[c$country=="US"] <- "US"
c$countryt[c$country=="CA"] <- "Canada*"
c$countryt[c$country=="JP"] <- "Japan"
c$countryt[c$country=="IS"] <- "Israel"
c$countryt[c$country=="SE"] <- "Senegal"
c$countryt[c$country=="DK"] <- "Denmark"
c$countryt[c$country=="IT"] <- "Italy*"
c$countryt[c$country=="UK"] <- "UK"
c$countryt[c$country=="FR"] <- "France*"
c$countryt[c$country=="IN"] <- "India"
c$countryt[c$country=="CN"] <- "China"
c$countryt[c$country=="NZ"] <- "New Zealand"
c$countryt[c$country=="SW"] <- "Sweden*"
c$countryt[c$country=="IN"] <- "India"
c$countryt[c$country=="CN"] <- "China"
myplot <- ggplot(c, aes(y=countryt, x=negcoef.rmssd)) +
  geom_density_ridges_gradient(rel_min_height=.05, scale=1.6, color="white", fill="black") +
  scale_fill_gradient(low="red",high="blue") + 
  xlab("The Impact of Negative News on RMSSD \n(Distributions of Individual-Level Coefficients, by Country)") + 
  ylab("") +
  scale_x_continuous(limits=c(-40,40)) +
  scale_y_discrete(limits=c("US","UK","Sweden*","Senegal","Russia","New Zealand","Japan","Italy*","Israel","India",
                            "Ghana","France*","Denmark","China","Chile","Canada*","Brazil*")) +
  theme(
    axis.ticks.y=element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.y = element_text(color = "black", size = 11, vjust = -.2),  
    axis.text.x = element_text(color = "black", size = 10),  
    panel.background = element_rect(fill = "white",colour = "white",size = 0.5, linetype = "solid"),
    panel.grid.major = element_line(size = 0.25, linetype = 'solid',colour = "gray"), 
    panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "white")
  ) +
  geom_vline(xintercept = 0, color="white") 
suppressWarnings(print(myplot))

#####
#GSL results in Figure 4

c <- sv[sv$dropgsl.new==0,c("negcoef.gsl","country")]
c <- na.omit(c)
c$countryt <- NA
c$countryt[c$country=="GH"] <- "Ghana*"
c$countryt[c$country=="RU"] <- "Russia*"
c$countryt[c$country=="BR"] <- "Brazil*"
c$countryt[c$country=="CH"] <- "Chile*"
c$countryt[c$country=="US"] <- "US*"
c$countryt[c$country=="CA"] <- "Canada*"
c$countryt[c$country=="JP"] <- "Japan"
c$countryt[c$country=="IS"] <- "Israel*"
c$countryt[c$country=="SE"] <- "Senegal*"
c$countryt[c$country=="DK"] <- "Denmark*"
c$countryt[c$country=="IT"] <- "Italy"
c$countryt[c$country=="UK"] <- "UK"
c$countryt[c$country=="FR"] <- "France"
c$countryt[c$country=="IN"] <- "India"
c$countryt[c$country=="CN"] <- "China"
c$countryt[c$country=="NZ"] <- "New Zealand"
c$countryt[c$country=="SW"] <- "Sweden"
c$countryt[c$country=="IN"] <- "India"
c$countryt[c$country=="CN"] <- "China"
myplot <- ggplot(c, aes(y=countryt, x=negcoef.gsl)) +
  geom_density_ridges(rel_min_height=.05, scale=1.6, color="white", fill="black") +
  xlab("The Impact of Negative News on nSCL \n(Distributions of Individual-Level Coefficients, by Country)") + 
  ylab("") +
  scale_x_continuous(limits=c(-.02,.02)) +
  scale_y_discrete(limits=c("US*","UK","Sweden","Senegal*","Russia*","New Zealand","Japan","Italy","Israel*","India",
                            "Ghana*","France","Denmark*","China","Chile*","Canada*","Brazil*")) +
  theme(
    axis.ticks.y=element_blank(),
    panel.grid.major.x = element_blank(),
    axis.text.y = element_text(color = "black", size = 11, vjust = -.2),  
    axis.text.x = element_text(color = "black", size = 10),  
    panel.background = element_rect(fill = "white",colour = "white",size = 0.5, linetype = "solid"),
    panel.grid.major = element_line(size = 0.25, linetype = 'solid',colour = "gray"), 
    panel.grid.minor = element_line(size = 0.25, linetype = 'solid',colour = "white")
  ) +
  geom_vline(xintercept = 0, color="white") 
suppressWarnings(print(myplot))

#####
#HRV and GSL results in Table S5

stata_src <- "
  xtset resp order
  quietly xtreg rmssd negativity order subtitles , re 
  estimates store A
  quietly xtreg rmssd c.negativity##i.subtitles order , re 
  estimates store B
  esttab A B, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w r2_b r2_o N N_g)
  "
stata(stata_src,data.in=rs.hrv)

stata_src <- "
  xtset resp order
  quietly xtreg gslmc negativity order subtitles, re 
  estimates store A
  quietly xtreg gslmc c.negativity##i.subtitles order , re 
  estimates store B
  esttab A B, cells(b(star fmt(3)) se(par fmt(3))) nogaps stat(r2_w r2_b r2_o N N_g)
  "
stata(stata_src,data.in=rs.gsl)


#####
#Table S6

cor.test(sv$negcoef.rmssd,sv$media,use="complete.obs")
cor.test(sv$negcoef.rmssd,sv$gdppc,use="complete.obs")
cor.test(sv$negcoef.rmssd,sv$hofstede,use="complete.obs")
cor.test(sv$negcoef.gsl,sv$media,use="complete.obs")
cor.test(sv$negcoef.gsl,sv$gdppc,use="complete.obs")
cor.test(sv$negcoef.gsl,sv$hofstede,use="complete.obs")

